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Abstract 

The standard hybrid Monte Carlo algorithm uses the second order integrator at 
the molecular dynamics step. This choice of the integrator is not always the best. 
Using the Wilson fermion action, we study the performance of the hybrid Monte 
Carlo algorithm for lattice QCD with higher order integrators in both zero and finite 
temperature phases and find that in the finite temperature phase the performance of 
the algorithm can be raised by use of the 4th order integrator. 
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1 Introduction 



In lattice QCD the hybrid Monte Carlo (HMC) algorithm |Q] is widely used for 
simulations of even- flavor of quarksQ These simulations are usually difficult tasks, 
especially at small quark masses where the computational cost of the matrix solver 
which is the most time consuming part of the HMC algorithm grows. In order to 
obtain reliable results within limited computer resources it is important to find an 
efficient way to implement the HMC algorithm so that the total computational cost 
of the algorithm is minimized || . 

The basic idea of the HMC is a combination of 1) molecular dynamics and 2) 
Metropolis test. Usually the second order leapfrog integrator is used at the molecular 
dynamics (MD) step. The integrator causes 0(At 3 ) integration errors, where At 
denotes the step-size. Due to the integration errors the Hamiltonian is not conserved. 
The errors introduced by this integrator have to be removed by the Metropolis test, 
i.e. accept the new configuration with a probability ~ min(l, exp(— AH)) where 
AH = H new — H id is an energy difference between the starting Hamiltonian H id 
and the new Hamiltonian H new at the end of the trajectory. 

The acceptance at the Metropolis step depends on the magnitude of the energy 
difference induced by the integration errors. If a higher order integrator is used at the 
MD step, the integration errors can be reduced. Therefore one may easily imagine 
that the performance of the HMC increases with the higher order integrator. However 
this is not always true since the higher order integrator has more arithmetic operations 
than the lower one and this might decrease the performance. The total performance 
should be measured by taking account of two effects: acceptance and the number of 
arithmetic operations. In Ref.[Q] the performance of the higher order integrator at 
zero temperature (/3 = 0) was studied systematically and it turned out that for the 
simulation parameters used for the current large-scale simulations with the Wilson 
fermion action the 2nd order integrator is the best one. The main reason why the 
higher order integrators are not so effective is that the energy difference caused by 
the errors of the higher order integrator increases more rapidly than that of the lower 
one as the quark mass decreases. In Ref.[||] it is shown that at finite temperature 
the quark mass dependence of the energy difference is small. If so, the conclusion 
of Ref.Q might change at finite temperature. In this letter we test higher order 
integrators at finite temperature and demonstrate that they can actually perform 
better than the lower order. 



Odd-flavor simulations of the HMC algorithm are also possible by modifying the Hamiltonian 
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2 Higher order integrator 




In this section we define higher order integrators studied here. Our definition is 
parallel to that of Ref.[[|]. Let H be a Hamiltonian given by 

H = \p 2 + S(q) (1) 

where q = (fc, ...) and p = (pi,P2,---) are coordinate variables and conjugate 
momenta respectively. S(q) is a potential term of the system considered. For the 
lattice QCD, q correspond to link variables and S(q) consists of gauge and fermion 
actions. 

In the MD step we solve Hamilton's equations of motion, 

_ dH 

dqi ' 

In general these equations are not solvable analytically. Introducing a step-size At, the 
discretized version of the equations are solved. In the conventional HMC simulations 
the 2nd order leapfrog scheme, which causes 0(At 3 ) step-size error, is used to solve 
the equations. This scheme is written as 

«(« + ¥) =?(t) + fp(*) 

p(t + At) = p(t )-A t ££Mi±f)) (3, 

oq 

q(t + At) = q{t + ^) + ^p{t + At). 

Eq.([|) forms an elementary MD step. This elementary MD step is performed repeat- 
edly iV times. The trajectory length r is given by r = N x At. 

Any integrator which satisfies two conditions: l)time reversible and 2)area pre- 
serving can be used for the MD step of the HMC. These conditions are needed to 
satisfy the detailed balance. When we use the Lie algebraic formalism || fj], ||, ||] we 
can easily construct higher order integrators which satisfy the above conditions. From 
the Lie algebraic formalism we find that higher order integrators can be constructed 
by combining lower order integrators. Let G2 n d(At) be an elementary MD step of 
the 2nd order integrator with a step-size At. The 4th order integrator Guh{At) is 
constructed from a product of three 2nd order integrators as ^, ||, ||, [l0| 

G 4th (At) = G 2 nd(aiAt)G2nd(a 2 At)G2nd(aiAt), (4) 

where the coefficients a« are given by 

2 l/3 

a2 = ~J^2^- (6) 
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Eq.([|) means that there are three elementary MD steps: (i) first we integrate the 
equations by eq.(||) with a step-size of a\At, (ii) then proceed the 2nd order integra- 
tion with a step-size of a 2 At, (iii) finally integrate the equations with a step-size of 
a\At. After performing these three elementary MD steps sequentially we obtain the 
4th order integrator with the step-size At. This construction scheme can be general- 
ized to the higher even-order integrators. The (2k + 2)-th order integrator is given 
recursively by 

G 2k+2 {At) = G 2k (b 1 At)G 2k (b 2 At)G 2k (b 1 At), (7) 
where the coefficient bi are 

bl = 2 - 2V(2*+D (8) 

2 l/(2fc+l) 

62 = ~ 2 -2 1 /(2fc+i)- ( 9 ) 

Compared with the 2nd order integrator, the computational cost of the n-th order 
one constructed from eq.(0) grows as 3 n//2_1 . For instance the 6th order integrator 
has 9 times more arithmetic operations than those of the 2nd order one. Yoshida || 
found three parameter sets of 6th order integrators with less arithmetic operations 
(Table.^j). Yoshida's 6th order integrators consist of seven G 2nd s as 

G 6th (At) = G 2nd (w 3 At)G 2nd (w 2 At)G 2nd ( Wl At)G 2nd (w At) 

xG 2nd ( Wl At)G 2nd (w 2 At)G 2nd (w 3 At). (10) 

In Ref.[|| these 6th order integrators were examined and one of three parameter sets, 
denoted by Yl, was found to give smaller integration errors than those of others. In 
this study the parameter set Yl ( see Table.[l] ) is used for the 6th order integrator. 





Yl 


Y2 


Y3 




-0.117767998417887e-l 


-0.213228522200144e+l 


0.152886228424922e-2 


w 2 


0.235573213359357e+0 


0.426068187079180e-2 


-0.214403531630539e+l 


w 3 


0.784513610477560e+0 


0.143984816797678e+l 


0.144778256239930e+l 



Table 1: Parameter sets (Y1-Y3) of the 6th order integrators by Yoshida ||]. wq is given 
Wq = 1 — w\ — w 2 — W3 . 



3 Efficiency of the HMC algorithm 

In order to compare among various integrators one needs a criterion which ranks 
integrators. Following Ref.Q we utilize the efficiency function Ejf constructed from 
a product of acceptance P acc and step-size At: 

E ff = P acc At. (11) 
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This function has one maximum at a certain step-size which we denote At op t- Using 
the energy difference AH the acceptance is expressed by || 

(P acc ) = erfc^Atf) 1 / 2 ), (12) 



where erfc is the complementary error function. In stead of using eq.(12), when (AH) 
is small, we may use 

(P acc )=exp(-A ( l Aj ff 2 ) 1 /2). (13) 

Although mathematically speaking eq.(13) is valid only for {AH 2 ) 1 / 2 <C 1, the 
numerical study @ shows that eq.(13) approximates the acceptance quit well for 
(AH 2 ) 1 / 2 < 3 which corresponds to (P acc ) > 20%. This is enough for our purpose 
since typically the acceptance of the HMC is taken to be (P aC c) > 50%. 

In the lowest order of At, (AH 2 ) 1 / 2 of the n-th order integrator is expressed as 

SI ' , 

(AH 2 ) 1 / 2 = C n V 1/2 At n + 0(At n+1 ), (14) 

where V is volume of the system and C n is a Hamiltonian dependent coefficient. 
Substituting eq.(|i~4|) into eq.(^) one obtains 

(P acc ) =eM~CnV^At n ), (15) 



where C n = C n /v 2vr. If one uses eq.(|15|) for eq.(|ll|) one can easily obtain the optimal 
step-size which gives the maximum efficiency: 



At opt = H— r- (16) 

V nC n V2 

Furthermore substituting At opt to eq. (|l5|) one obtains the optimal acceptance^] as 0] 

(Pacc)opt = exp(— ) (17) 
n 

0.61 2nd 

0.78 4th (18) 
0.85 6th. 

Note that the above result does not depend on the specific Hamiltonian and can be 
applied for any model. 

Using eq.(|i"6l) and ( |l7l ) we obtain the optimal efficiency of the n-th order integrator: 



( E ff)ort th = eM-kd^-T, (19) 
n V nC n V2 



2 A recent study |Tl| shows that if one considers higher order effects of At the optimal acceptance might 
be slightly changed. In the current study we omit this small effect. 
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We use eq.(|19D to compare among integrators. Eq,fll9|) is easy to handle because 
eq. (|l9|) has one unknown parameter C n and the value of C n can be estimated easily 
from a single simulation on a rather small lattice. 

Now let us compare n-th and m-th order integrators (n > m). If one obtains a 
gain G for the n-th order integrator over the m-th order one, the following condition 
must be satisfied: 

( E ff)opt th = Gk nm (E ff )™t th , (20) 

where k nm is a relative cost factor needed to implement the n-th order integrator 
against the m-th one and G is defined so that both n-th and m-th order integrators 
are equally effective with G = 1. In our case, /C42 = 3 and k$4 = 7/3. Substituting 



eq.(19) to eq.(^) and rewriting the equation one obtains the lattice volume needed 



to have the gain G: 

(\ n 
— L (nC n ) m . (21) 
mC m J 

Our formulas are based on the assumption that (AH 2 ) 1 / 2 satisfies eq.(14). The 



validity of eq.(14) depends on the action we take and simulation parameters. In 
general it is expected that the contribution of the higher order terms to eq.(14) 
becomes small as the lattice size increases. Let us rewrite eq,(|l4|) as 

(AH 2 ) 1 / 2 = C n V l / 2 At n + C n+l V 1 / 2 At n+1 , (22) 

where Cn+iV 1 / 2 At n+1 stands for the first relevant higher order term. We measure 
the contribution of the higher order term by the ratio: 

Cn+lV^Af^ _ C n+1 

c n yi/ 2 At" - ~c7 AL {2d) 

In order to keep the constant acceptance, At must be taken to be small as the lattice 
size increases. Thus eq. (p3[) , ie. the contribution of the higher order term, becomes 
small for larger lattices. 

Fig.l shows a typical example of C n on a 16 3 x 4 lattice as a function of At. 
{AH 2 ) 1 / 2 are well expressed by the functions proportional to At™ except for large 



At. We are only interested in the acceptance region indicated by eq.(18): P acc > 
which corresponds to {AH 2 ) 1 / 2 < l|. It seems that for {AH 2 ) 1 / 2 < 1 the effect of 
the higher order terms is small. 

4 Simulation results 

We use the plaquette gauge action and standard Wilson fermion action with two 
flavors of degenerate quarks (n/ = 2). We first determine the coefficients C n at zero 

3 See Fig.l in f§. 
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and finite temperature. This can be done by measuring (AH 2 ) 1 / 2 at a small step-size 
and substituting the results into eq.(li~4|). The trajectory length r is set to 1. We 
choose (3 = 5.75 and make simulations on both 12 4 and 12 3 x 4 lattices. The critical 
kappa k c of rif = 2 at this (3 is around 0.157Q. We use several k's in a range of 
0.1 < k < 0.155. In this range we maintain the confinement phase on the 12 4 lattice 
and the deconfinement phase on the 12 3 x 4 lattice. In this study we refer to the 
results on the 12 4 (12 3 x 4) lattice as those at zero(finite) temperature. 

Figs. 2-4 show C n as a function of m q a. Here m q a is defined by m q a = (1/k — 
l/« c )/2. At large quark masses where the fermionic effects are negligible, values of 
C n at zero and finite temperature coincide each other. This may indicate that the 
contribution to C n from the gauge sector is almost same at zero and finite temper- 
ature. At small quark zero temperature increases more rapidly than 
those at finite temperature as m q a decreases. The quark mass dependence of C n at 
finite temperature, compared to that at zero temperature, is small for all the integra- 
tors studied here. This behavior is consistent with the result of Ref.[||. Substituting 
values of C n at finite temperature into eq.(^ll) with G = 1, we calculate the lattice 
size L e ( here V = L 3 x Nt )• This lattice size L e is the one with which the higher 
order integrator and the lower order one perform equally. For a lattice size L > L e 
the higher order integrator is more effective than the lower one. Fig. 5 shows L e from 
comparison between the 2nd and 4th order integrators (2nd vs. 4th) and between 
the 4th and 6th order ones (4th vs. 6th). For the case of 4th vs. 6th, L e increases 
as m q a decreases, which we do not appreciate. On the other hand for the 2nd vs. 
4th, L e remains less than 20 even at small m q a. This result is contrast to that ob- 
tained at zero temperature where L e increases as m q a decreases @. The above result 
encourages us to use the 4th order integrator at finite temperature. 

The results in Fig. 5, however, just show the lattice on which the both integrators 
are equally effective. To use the higher order integrator in simulations one must 
obtain some gain over the lower order one. In Fig. 6, using eq.(^), we show the 
expected gain ( the region between the solid lines ) at k = 0.1525(m g a ~ 0.094) as a 
function of lattice size L. To have G = 2 gain ( which means 2 times faster ) a lattice 
size L ~ 100 is required. This huge lattice size is still not accessible in the current 
large-scale simulations. Probably the maximum lattice size accessible at the moment 
is L < 50. Therefore we can not expect a large gain from the 4th order integrator 
even if it is used now. If we use a lattice with L ~ 50, G ~ 1.5 can be achieved. Thus 
at the level of the current large-scale simulations, we expect to obtain G < 1.5. 

We also make simulations at k = 0.1525 to confirm that we can actually obtain 
some gain for the 4th order integrator over the 2nd one. At k = 0.1525, L e is 
estimated to be 17.5(9). We choose 18 3 x 4 and 28 3 x 4 lattices. On the 18 3 x 4 lattice 
we expect G ~ 1 and on the 28 3 x 4 lattice, G > 1. The step-size is adjusted so that 

4 This value is estimated from Figure 5 in Ref. |l2| 
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18 3 x 4 


28 3 


x 4 


k = 0.1525 


2nd 4th 


2nd 


4th 


At 


1/24 1/10 


1/36 


1/12 


Acceptance 


0.57(2) 0.81(2) 


0.60(3) 


0.81(3) 



Table 2: Step-size and acceptance at k = 0.1525. 



the acceptance gives a similar value with eq.([TJ). The gain G is calculated by 

q = (Pace X At) At h ^ 
3(P a cc X At)2nd 

where a factor of 3 in the denominator comes from the relative cost factor = 3. 
Table 2 shows the simulation results. Using these results, we obtain G = 1.08(3) on 
the 18 3 x 4 lattice and G = 1.35(8) on the 28 3 x 4 lattice (see also Fig.6). As expected 
the gain increases with L. The result on the 28 3 x 4 lattice is an example showing 
that the 4th order integrator is more effective than the 2nd order one. 



5 Conclusions 

We have studied higher order integrators for the HMC algorithm with the Wilson 
fermion action at finite temperature. Contrast to the zero temperature case, the 4th 
order integrator at finite temperature can be more effective than the 2nd order one. 
This was demonstrated by the simulations at (3 = 5.75 on L 3 x 4 lattices. The gain 
is dependent of the lattice size. It was shown that on the 28 3 x 4 lattice at (3 = 5.75 
and k = 0.1525 the 4th order integrator is about 35% faster than the 2nd order one. 

When large-scale simulations at finite temperature are planned, it is recommended 
to check which integrator is effective for the lattice considered. This check can be done 
easily. First we measure C2 and C4. This first step does not take much computational 
time since they can be measured on a small lattice. Then substituting values of C2 
and C4 to eq.(|2l|), we obtain a relation between G and L. If we obtain G > 1 on the 
lattice planned for the simulations, we should consider to use the 4th order integrator. 

Eq.(^T|) is obtained by using the approximation at small At. There might exist 
the systematic errors caused by the approximation. Although for the present study 
we considered the errors to be small, for other actions and simulation parameters 
they might contribute largely. We have to keep in mind that there might exist such 
contributions depending on the simulation details. 
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Figure 1: (AH 2 ) 1 / 2 on a 16 3 x 4 lattice at j3 = 5.75 and k = 0.155 as a function of At. 
The lines proportional to At n are also drawn. The simulations were done with the Wilson 
fermion action. 
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Figure 2: C2 as a function of the quark mass m q a. In the figure, C2 is normalized as 
C 2 /(2vr) 1 /2 = c 2 . 
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Figure 3: Same as in Fig. 2 but for C4. 
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Figure 4: Same as in Fig. 2 but for Cq. 
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Figure 5: L e as a function of the quark mass m q a. L e stands for the lattice size for 
which the higher order integrator and the lower order one are equally effective. Circles are 
from comparison between the 2nd order integrator and the 4th order one, and squares are 
comparison between the 4th order integrator and the 6th order one. 
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Figure 6: The gain G as a function of L ( here V = L 3 x 4 ) at k = 0.1525. The solid 
lines are calculated from eq . fl2l"D . Circles are from Monte Carlo simulations. 
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